C
      SUBROUTINE GCPFIT (GCP, MGCP, MAXDEG)
C
C     TO FIND GEOMETRIC TRANSFORMATION NEEDED FOR GEOGRAPHIC REFERENCING
C
      DIMENSION GCP(4,MGCP)
      DOUBLE PRECISION COEF(21,4)
      LOGICAL CCT
      INTEGER NTERM(5)/3,6,10,15,21/
      COMMON /LANDST/ CCT, LLC, DEGCEN, SAMPOF, LINOFF, AMPL, PHASE
      COMMON /LSQCFC/ COEF, LSQDEG, IER1(2), TOL
      COMMON / MEANS/ GCPM(4)
C
C     GCP IS GROUND CONTROL POINT TABLE INPUT BY USER
C     1.  PIXEL WITHIN LINE
C     2.  LINE WITHIN FRAME
C     3.  EASTING
C     4.  NORTHING
C     COEF(21,I)- COEFFICIENTS FOR : 1- PIXEL
C                                    2- LINE
C                                    3- EASTING
C                                    4- NORTHING
C
      IF (CCT) CALL GCPCOR (GCP, MGCP)
      AMGCP=MGCP
      RADDEG = 180.0/3.14159265
      TOL=1.E-20
      DO 1 I=1,4
    1 GCPM(I) = 0.0
      WRITE (6,101)
      WRITE (6,1020)
C
C     COMPUTE SUMS OF INPUT DATA
      DO 10 J=1,MGCP
      DO 5 I=1,4
    5 GCPM(I) = GCPM(I) + GCP(I,J)
      WRITE (6,102) J, (GCP(I,J), I=1,4)
   10 CONTINUE
C
C     COMPUTE MEANS OF INPUT DATA
      DO 11 I=1,4
   11 GCPM(I) = GCPM(I) / AMGCP
      WRITE (6,104) GCPM
C
C     SUBTRACT MEANS OF INPUT DATA
      WRITE (6,1019)
      WRITE (6,1020)
      DO 13 J=1,MGCP
      DO 12 I=1,4
   12 GCP(I,J) = GCP(I,J) - GCPM(I)
      WRITE (6,102) J, (GCP(I,J), I=1,4)
   13 CONTINUE
C
C     FIND POLYNOMIAL FITS FOR DEGREES 1 - INPUT VALUE OF LSQDEG
      DO 100 LSQDEG=1,MAXDEG
      ISTOP=NTERM(LSQDEG)
      IF (ISTOP.GT.MGCP) GO TO 400
      CALL LSQCF (0, GCP, MGCP)
      WRITE(6,383) LSQDEG,IER1,((COEF(J,I),I=1,4),J=1,ISTOP)
      WRITE(6,111)
      RESPX=0.0
      RESUT=0.0
C
C     PERFORM ACCURACY ANALYSIS OF FIT
      DO 388 I=1,MGCP
      GP=GCP(1,I) + GCPM(1)
      GL=GCP(2,I) + GCPM(2)
      GE=GCP(3,I) + GCPM(3)
      GN=GCP(4,I) + GCPM(4)
C
      CALL EVPOLY(1,GE,GN,ANS)
      DP=GP-ANS
      CALL EVPOLY(2,GE,GN,ANS)
      DL=GL-ANS
      SQ = DP**2 + DL**2
      RESPX = RESPX + SQ
      XMAG = SQRT (SQ)
      XDIR = RADDEG * ATAN2(-DL,DP)
C
      CALL EVPOLY(3,GP,GL,ANS)
      DE=GE-ANS
      CALL EVPOLY(4,GP,GL,ANS)
      DN=GN-ANS
      SQ = DE**2 + DN**2
      RESUT = RESUT + SQ
      UMAG = SQRT (SQ)
      UDIR = RADDEG * ATAN2(DE,DN)
      WRITE (6,109) I, XMAG, XDIR, DP, DL, UMAG, UDIR, DE, DN
388   CONTINUE
C
      AGCP1=MGCP-1
      RESPX = SQRT (RESPX/AGCP1)
      RESUT = SQRT (RESUT/AGCP1)
      WRITE (6,390) RESPX, RESUT
  100 CONTINUE
      LSQDEG = MAXDEG
      RETURN
C
  400 WRITE (6,1100) LSQDEG, ISTOP, MGCP
      LSQDEG = LSQDEG - 1
      RETURN
C
C     FORMAT STATEMENTS.
  101 FORMAT ('1',20X,'GROUND CONTROL POINTS'/)
  102 FORMAT (1X,I4,4F15.3)
  104 FORMAT (///20X,'MEANS OF INPUT DATA ARE'//5X,4F15.3/
     .     '0THESE MEANS ARE FIRST SUBTRACTED')
  109 FORMAT (1X,I3,2(F10.3,F10.1,8X,2F10.3,8X))
  111 FORMAT (/35X,'COMPARISON OF OBSERVED AND PREDICTED VALUES'/
     .25X,'LANDSAT',50X,'GEOGRAPHIC'/15X,'ERROR',20X,'OBS - PRED',20X,
     .'ERROR',20X,'OBS - PRED'/5X,'MAGNITUDE  DIRECTION',10X,'P ERROR
     .L ERROR',10X,'MAGNITUDE  DIRECTION',10X,'E ERROR   N ERROR'/)
  383 FORMAT ('1  LEAST SQUARES FIT OF DEGREE',I3,25X,'ERROR CODES =',
     .2I3/38X,'COEFFICIENTS'/11X,'PIXEL',15X,'LINE',15X,'EASTING',12X,
     .'NORTHING'/(1X,1P4E20.6))
  390 FORMAT (/2(15X,'RMS ERROR =',E13.6,15X))
 1019 FORMAT ('1',7X'CONTROL POINTS AFTER SUBTRACTING MEANS'/)
 1020 FORMAT (20X,'LANDSAT',20X,'GEOGRAPHIC'/18X,'COORDINATES',18X,'COOR
     .DINATES'/' GCP NO.',7X,'PIXEL',10X,'LINE',9X,'EASTING',7X,'NORTHIN
     .G'/)
 1100 FORMAT (/20X,'FIT OF DEGREE',I2,' REQUIRES',I3,' GROUND CONTROL PO
     .INTS.',I5,' WERE SUPPLIED.'/)
      END
C
      SUBROUTINE LSQCF (IDEL, GCP, MGCP)
C
C     THIS SUBROUTINE CALCULATES THE LEAST-SQUARES COEFFICIENTS
C     FOR THE BIVARIATE POLYNOMIAL
C
      DIMENSION GCP(4,MGCP)
      DOUBLE PRECISION COEF(21,4), CORE(2500), X(45), B(200), AUX(100),
     .A(2100), S(45), XXX, YYY, AMAX, C
      INTEGER NTERM(5)/3,6,10,15,21/
      INTEGER XP(21)/0,1,0,2,1,0,3,2,1,0,4,3,2,1,0,5,4,3,2,1,0/
      INTEGER YP(21)/0,0,1,0,1,2,0,1,2,3,0,1,2,3,4,0,1,2,3,4,5/
      INTEGER IPIV(50),IND1(2)/1,3/,IND2(2)/3,1/
      EQUIVALENCE (CORE(1),X(1)), (CORE(46),B(1)), (CORE(246),AUX(1)),
     .(CORE(346),A(1)), (CORE(2446),S(1))
      COMMON /LSQCFC/ COEF, LSQDEG, IER1(2), TOL
C
C     GCP IS GROUND CONTROL POINT TABLE INPUT BY USER
C     1.  PIXEL WITHIN LINE
C     2.  LINE WITHIN FRAME
C     3.  EASTING
C     4.  NORTHING
C     COEF(21,I)- COEFFICIENTS FOR : 1- PIXEL
C                                    2- LINE
C                                    3- EASTING
C                                    4- NORTHING
      L=MGCP
      L1=MGCP
      IF(IDEL.NE.0) L=MGCP-1
      ISTOP=NTERM(LSQDEG)
      IF (ISTOP.GT.L) RETURN
C
      DO 10 II=1,2
      III=IND1(II)
      IIII=(II-1)*2+1
      IV=IND2(II)
      NGCP1=0
C
      DO 50 NGCP=1,L1
      IF(NGCP.EQ.IDEL) GO TO 50
      NGCP1=NGCP1+1
      B(NGCP1)=GCP(III,NGCP)
      NN=NGCP1+L
      B(NN)=GCP(III+1,NGCP)
      DO 50 I=1,ISTOP
      XXX=1.0D0
      IF(XP(I).NE.0) XXX=GCP(IV,NGCP)**XP(I)
      YYY=1.0D0
      IF(YP(I).NE.0) YYY=GCP(IV+1,NGCP)**YP(I)
      K=(I-1)*L+NGCP1
      A(K)=XXX*YYY
50    CONTINUE
C
C     SCALE THE MATRIX
      DO 70 I=1,ISTOP
      AMAX=0.D0
C
      DO 60 J=1,L
      K=(I-1)*L+J
      C=DABS(A(K))
60    AMAX=DMAX1(AMAX,C)
      IF(AMAX.EQ.0.D0) AMAX=1.D0
C
      DO 70 J=1,L
      K=(I-1)*L+J
      A(K)=A(K)/AMAX
70    S(I)=AMAX
C
      CALL DLLSQ(A,B,L,ISTOP,2,X,IPIV,TOL,IER,AUX)
      IER1(II)=IER
      DO 80 I=1,ISTOP
      J=I+ISTOP
      COEF(I,IIII)=X(I)/S(I)
      COEF(I,IIII+1)=X(J)/S(I)
80    CONTINUE
10    CONTINUE
      RETURN
      END
C
      SUBROUTINE EVPOLY (IFUN, X, Y, ANS)
C
C     EVALUATE POLYNOMIAL FIT FUNCTIONS
C
      DOUBLE PRECISION COEF(21,4), ANSD, XD, YD, XXX, YYY
      INTEGER NTERM(5)/3,6,10,15,21/
      INTEGER XP(21)/0,1,0,2,1,0,3,2,1,0,4,3,2,1,0,5,4,3,2,1,0/
      INTEGER YP(21)/0,0,1,0,1,2,0,1,2,3,0,1,2,3,4,0,1,2,3,4,5/
      COMMON /LSQCFC/ COEF, LSQDEG, IER1(2), TOL
      COMMON / MEANS/ GCPM(4)
C
C     COEF(21,I)- COEFFICIENTS FOR : 1- PIXEL
C                                    2- LINE
C                                    3- EASTING
C                                    4- NORTHING
      ISTOP = NTERM(LSQDEG)
      ANSD = 0.0
      IF (IFUN.EQ.3.OR.IFUN.EQ.4) GO TO 1
      XD = X - GCPM(3)
      YD = Y - GCPM(4)
      GO TO 5
    1 XD = X - GCPM(1)
      YD = Y - GCPM(2)
    5 CONTINUE
C
      DO 10 I=1,ISTOP
      XXX = 1.0
      IF (XP(I).NE.0) XXX = XD**XP(I)
      YYY = 1.0
      IF (YP(I).NE.0) YYY = YD**YP(I)
   10 ANSD = ANSD + COEF(I,IFUN)*XXX*YYY
C
      ANS = ANSD + GCPM(IFUN)
      RETURN
      END
C
      SUBROUTINE GCPCOR (GCP, MGCP)
C
C     GCPCOR REMOVES MIRROR VELOCITY PROFILE AND EARTH CURVATURE
C     DISTORTIONS PRIOR TO OBTAINING LEAST SQUARES FITS
C
C     GCP IS GROUND CONTROL POINT TABLE INPUT BY USER
C     1.  PIXEL WITHIN LINE
C     2.  LINE WITHIN FRAME
C     3.  EASTING
C     4.  NORTHING
C
      DIMENSION GCP(4,MGCP)
      REAL*4 MVPOFF
      COMMON /LANDST/ CCT, LLC, DEGCEN, SAMPOF, LINOFF, AMPL, PHASE
C
C     ADD BACK DELAY TO GET ORIGINAL CCT PIXEL NUMBER
      WRITE (6,100)
      DO 10 NGCP = 1,MGCP
      WRITE (6,101) GCP(2,NGCP), GCP(1,NGCP)
      LINE = GCP(2,NGCP) + 0.5
      SAMP = GCP(1,NGCP) + DELAY(LINE)
C
C     APPLY MIRROR VELOCITY PROFILE CORRECTION
      OFF = MVPOFF(SAMP)
      GCP(1,NGCP) = GCP(1,NGCP) - OFF
C
C     APPLY EARTH CURVATURE CORRECTION
      CURV = ERCURV(SAMP-OFF)
      GCP(1,NGCP) = GCP(1,NGCP) - CURV
C
      SAMP0 = SAMP + SAMPOF
      SHIFT = OFF + CURV
      WRITE (6,102) SAMP, SAMP0, OFF, CURV, SHIFT, GCP(1,NGCP)
   10 CONTINUE
      RETURN
C
  100 FORMAT ('1',20X,'GROUND CONTROL POINT CORRECTIONS'//19X,'OFFSET',
     .5X,'ORIGINAL',3X,'CCT SCENE',7X,'MVP',9X,'CURV',8X,'NET',8X,
     .'CORR.'/7X,'RECORD',3(6X,'SAMPLE'),3(7X,'SHIFT'),6X,'OFFSET'/)
  101 FORMAT (1X,2F12.3)
  102 FORMAT ('+',24X,8F12.3)
      END
C
      FUNCTION DELAY (LINE)
C
C     COMPUTE ROTATIONAL AND SENSOR DELAY RELATIVE TO FIRST SWATH
C
C     DEGCEN - LATITUDE AT THE CENTER OF THE LANDSAT SCENE
C     PIXDLY - EQUATORIAL EARTH ROTATION PER SWATH IN PIXELS
C     DEGPLN - CHANGE IN LATITUDE PER SCAN LINE
C     RADDEG - RADIANS PER DEGREE
C     SDELAY - SENSOR SAMPLING INTERVAL BETWEEN LINES (2) / TOTAL NUMBER
C              OF DETECTORS (25)
C     LINESW - LINE NUMBER IN THE SWATH (1 - 6)
C
      DIMENSION SDELAY(6)
      INTEGER CCTLIN, SWATH
      LOGICAL CCT
      COMMON /LANDST/ CCT, LLC, DEGCEN, SAMPOF, LINOFF, AMPL, PHASE
      DATA PIXDLY, DEGPLN, RADDEG, SDELAY /0.6, 0.00071086, 0.01745329,
     .0.0, 0.08, 0.16, 0.24, 0.32, 0.40/
C
      IF (CCT) GO TO 10
      DELAY=0.0
      RETURN
C
   10 CONTINUE
      CCTLIN = LINE + LINOFF
      DEGLAT = DEGCEN + (1170.5-CCTLIN) * DEGPLN
      RDELAY = PIXDLY * COS (DEGLAT*RADDEG)
      SWATH = (CCTLIN-1)/6 + 1
      LINESW = MOD(CCTLIN-1,6) + 1
      DELAY = RDELAY*(SWATH-1) - SDELAY(LINESW)
      RETURN
      END
C
      REAL FUNCTION MVPOFF (PS)
C
C     COMPUTE MIRROR VELOCITY PROFILE OFFSET
C
C     LLC - NUMBER OF PIXELS IN THE RAW SCAN LINE
C     SAMPOF - NUMBER OF PIXELS SKIPPED IN THE LANDSAT SCENE
C     AMPL, PHASE - AMPLITUDE, PHASE OF MIRROR VELOCITY PROFILE CURVE
C
      LOGICAL CCT
      COMMON /LANDST/ CCT, LLC, DEGCEN, SAMPOF, LINOFF, AMPL, PHASE
C
      IF (CCT) GO TO 12
      MVPOFF = 0.0
      RETURN
C
   12 CONTINUE
      PS1 = PS + SAMPOF
      MVPOFF = AMPL * SIN (6.2831853 * (PS1+PHASE-1.0) / (LLC-1))
      RETURN
      END
C
      FUNCTION ERCURV (PS)
C
C     COMPUTE EARTH CURVATURE CORRECTION
C
C     RE, RSAT - EARTH RADIUS, SATELLITE ORBIT RADIUS
C     TFOV - TOTAL FIELD OF VIEW OF THE SCANNER (11.56 DEGREES)
C     LLC  - NUMBER OF PIXELS IN THE RAW SCAN LINE
C     SAMPOF - NUMBER OF PIXELS SKIPPED IN THE LANDSAT SCENE
C
      LOGICAL CCT
      COMMON /LANDST/ CCT, LLC, DEGCEN, SAMPOF, LINOFF, AMPL, PHASE
      DATA RE, RSAT, TFOV /6367.4, 7285.6, 0.20176/
C
      IF (CCT) GO TO 10
      ERCURV = 0.0
      RETURN
C
C     COMPUTE SCANNER ANGLE AT PIXEL NO. PS AND ANGLE SUBTENDED AT
C     THE CENTER OF THE EARTH
   10 CONTINUE
      TFOV2 = TFOV/2.0
      PS1 = PS + SAMPOF
      ANGSCN = (PS1-1.0)*TFOV/(LLC-1) - TFOV2
      ANGERT = ARSIN(SIN(ANGSCN)*RSAT/RE) - ANGSCN
      TOTERT = ARSIN (SIN(TFOV2)*RSAT/RE) - TFOV2
C
C     FIND SCANNER ANGLE BASED ON FRACTION OF TOTAL ARC ON THE EARTH'S
C     SURFACE AND CONVERT TO PIXEL NUMBER
      ANGSC1 = TFOV2 * ANGERT / TOTERT
      PS2 = 1.0 + (ANGSC1+TFOV2) * (LLC-1) / TFOV
      ERCURV = PS1 - PS2
      RETURN
      END
